We simulated a trait controlled by 1,000 QTL 100 times. For each, one population was selected for an increase in the trait and one population was not selected at all (drift only). Selection or drift persisted for a variable numebr of generations, with a variable number of post-selection generations simulated. Herein, we use these simulations to test the power of our test statistic to identify selection on the trait.
These values will store results of the loop
Pvals.20.5 <- c()
source("Ghat.R")
Simulations were performed on the IRCF Biocluster using the software QMSim. Each iteration of the loop will load data, manipulate data, conduct a test for selection, and record the p-value.
for(sim in 1:100){
iter <- sprintf("%03i",sim)
directory="/home/beissingert/data/PolySelect_bigsims/Intermittent/r_intermittent"
#load genotypes and phenotypes
pheno <- read.table(paste(directory, "/20Select_5Not_data_",iter,".txt",sep=""),header=T,stringsAsFactors=F)
map <- read.table(paste(directory, "/lm_mrk_",iter,".txt",sep=""),header=T,stringsAsFactors=F)
geno <- read.table(paste(directory, "/20Select_5Not_mrk_",iter,".txt-head",sep=""),header=F,stringsAsFactors=F,skip=1,sep="", colClasses=c("numeric","character")) # read genos as characters #
qtlMap <- read.table(paste(directory, "/lm_qtl_",iter,".txt",sep=""),header=T,stringsAsFactors=F)
#load allele frequencies
freqs0<-read.table(paste(directory, "/20GenSelect_freq_mrk_",iter,".txt",sep=""),header=T,nrows=nrow(map),fill=T,stringsAsFactors=F) #NOTE POP TAKEN FROM
freqs20<-read.table(paste(directory, "/20Select_5Not_freq_mrk_",iter,".txt",sep=""),header=F,skip={nrow(map)+1},fill=T,stringsAsFactors=F) #
### Manipulate genotypes by coding to -1,0,1 and so that markers are columns, individuals are rows.
gen <- matrix(NA,nrow=nrow(map),ncol=nrow(geno))
gen <- as.data.frame(gen)
names(gen) <- geno[,1]
for(i in 1:1000){
#print(i)
tmp <- as.numeric(unlist(strsplit(geno[i,2],split="")))
tmp[which(tmp == 0)] <- -1
tmp[which(tmp == 3 | tmp ==4)] <- 0
tmp[which(tmp==2)] <- 1
gen[,i] <- tmp
}
gen<-t(gen)
gc()
### Calculate allele frequencies
#Gen 0
names(freqs0)[4]<- "Allele1"
names(freqs0)[5]<- "Allele2"
freqs0$Allele2[which(substr(freqs0$Allele1,1,1)==2)] <- "2:1.00000" # put this in the spot for the second allele
freqs0$Allele1[which(substr(freqs0$Allele1,1,1)==2)] <- "1:0.00000"
freqs0$Allele2[which(substr(freqs0$Allele1,3,3)==1)] <- "2:0.00000" ##
freqs0$Allele1 <- as.numeric(substr(freqs0$Allele1,3,1000))
freqs0$Allele2 <- as.numeric(substr(freqs0$Allele2,3,1000))
#Gen 20
names(freqs20)[4]<- "Allele1"
names(freqs20)[5]<- "Allele2"
freqs20$Allele2[which(substr(freqs20$Allele1,1,1)==2)] <- "2:1.00000" # put this in the spot for the second allele
freqs20$Allele1[which(substr(freqs20$Allele1,1,1)==2)] <- "1:0.00000"
freqs20$Allele2[which(substr(freqs20$Allele1,3,3)==1)] <- "2:0.00000" ##
freqs20$Allele1 <- as.numeric(substr(freqs20$Allele1,3,1000))
freqs20$Allele2 <- as.numeric(substr(freqs20$Allele2,3,1000))
#Calculate change
change2<-freqs20$Allele2-freqs0$Allele2
### Setup for function
geno <- gen
phen <- as.matrix(pheno[1:1000,10])
rownames(phen)<-pheno[1:1000,1]
change=change2
perms <- 10000
blockSize=100000/100
### Determine num_eff
potentials <- c(40000,13333,6667,4000,2857,2222,1818,1538,1333,1176,1053,667,400,286,222)
ldTab <- read.table(paste(directory,"/20Select_5Not_ld_decay_",iter,".txt",sep=""),skip=43,nrows=15)
ld <- as.numeric(substr(ldTab[,2],1,6))
Meff <- potentials[which(ld <= 0.03)[1]]
if(is.na(Meff)) Meff <- potentials[15] #
print(Meff)
# Run function
test<- Ghat_func(geno=geno,phen=phen,change=change2,method = "scale", num_eff = Meff, perms=1000,plot="Both")
Pvals.20.5[as.numeric(iter)] <- test$p.val
}
## [1] 400
## Loading required package: rrBLUP
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 667
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 667
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 667
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 286
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
## [1] 400
## [1] "Doing Permutation Test"
Below are some summary statistics about the results. Note also that I save the workspace.
hist(Pvals.20.5)
print(Pvals.20.5)
## [1] 1.382589e-05 2.337717e-08 1.131123e-05 5.328571e-09 6.049600e-08
## [6] 4.500957e-08 1.477593e-04 1.534103e-08 8.376485e-06 2.479922e-06
## [11] 6.809595e-07 4.008208e-08 7.051275e-06 4.410218e-09 1.323629e-08
## [16] 1.623001e-05 1.186751e-09 2.247514e-07 1.995555e-05 3.057348e-07
## [21] 1.385703e-05 7.575213e-09 5.230086e-06 1.994869e-10 7.539975e-06
## [26] 6.406634e-08 4.089458e-09 2.423449e-09 7.493862e-08 3.641873e-12
## [31] 3.225783e-05 1.894860e-05 2.882265e-07 1.828686e-08 7.477186e-05
## [36] 1.606267e-05 1.640856e-07 1.829604e-06 5.793842e-07 2.244819e-05
## [41] 7.061032e-06 9.341524e-06 2.058416e-05 2.861684e-04 6.405364e-05
## [46] 6.664228e-07 2.249467e-07 2.094596e-06 2.370846e-09 5.173802e-09
## [51] 2.278162e-06 2.602849e-06 3.985291e-10 9.442805e-07 1.230474e-05
## [56] 1.068580e-08 4.882379e-07 1.600186e-04 4.254955e-07 2.918771e-09
## [61] 6.553425e-06 1.262617e-06 1.368520e-07 2.047641e-06 4.820260e-06
## [66] 2.257799e-05 1.556989e-04 4.156825e-12 3.087850e-04 1.422297e-03
## [71] 5.909762e-05 7.818795e-06 1.956025e-07 4.274564e-07 2.178949e-06
## [76] 1.037955e-07 2.435861e-05 7.090692e-05 5.432631e-06 5.187017e-05
## [81] 3.667993e-04 4.589323e-05 7.803839e-07 3.726315e-08 1.658026e-06
## [86] 1.833542e-06 2.319005e-06 8.725873e-04 2.980138e-05 4.504325e-07
## [91] 4.490898e-06 8.184397e-05 2.504850e-03 6.434933e-07 3.618112e-09
## [96] 2.089820e-08 8.721400e-09 4.695622e-05 1.039105e-07 2.740904e-03
length(which(Pvals.20.5<=0.05))/length(Pvals.20.5)
## [1] 1
save.image("on20off5.RData")